home *** CD-ROM | disk | FTP | other *** search
/ TeX 1995 July / TeX CD-ROM July 1995 (Disc 1)(Walnut Creek)(1995).ISO / graphics / gnuplot / contrib / vanzandt / spline.c < prev    next >
C/C++ Source or Header  |  1992-10-25  |  12KB  |  444 lines
  1. /*    spline - interpolate points in a tabulated function or curve
  2.  
  3.     Usage...
  4.         spline [file][options]
  5.         options are:
  6.           -a  [step [start]] automatic abscissas 
  7.           -b                 break interpolation at each label
  8.           -c                 general curve 
  9.           -n  num            interpolate over num intervals 
  10.           -t  num            tension in interpolating curve 
  11.           -x  [min [max]]    interpolate from min to max 
  12.  
  13.     Notes...
  14.         This program fits a smooth curve through a given set of points,
  15.         using the splines under tension introduced by Schweikert [J. 
  16.         Math.  and Physics v 45 pp 312-317 (1966)].  They are similar
  17.         to cubic splines, but are less likely to introduce spurious
  18.         inflection points.
  19.  
  20.     History...
  21.         21 Sep 85    Added -xl and -yl switches for taking logs
  22.         23 Nov 85    Passing lines starting with ';' unchanged, otherwise
  23.                     ignoring them.
  24.  
  25.     Author...
  26.         Copyright (c) 1985  James R. Van Zandt
  27.  
  28.         All rights reserved.
  29.         This program may be copied for personal, non-profit use only.
  30.  
  31.         Based on algorithms by A.  K.  Cline [Communications of the ACM
  32.         v 17 n 4 pp 218-223 (Apr 74)].
  33.  
  34. */
  35.  
  36. #include <stdio.h>
  37. #include <math.h>
  38.  
  39. #define VERSION "1.2"
  40.  
  41. #define ENTRIES 200
  42. #define MAXLABELS 50
  43.  
  44. double *x, *y, *temp, *xp, *yp, *path;
  45. double slp1=0.,slpn=0.,sigma=-1.;
  46. double abscissa=0.;
  47. double abscissa_step=1.;
  48. double xmin=0.;
  49. double xmax=0.;
  50. double curv2(), mylog();
  51.  
  52. int xlog=0;        /* nonzero if taking log of x values */
  53. int ylog=0;        /* nonzero if taking log of y values */
  54. int debugging=0;
  55. int breaking=0;        /* nonzero if breaking at labels */
  56. int labels=0;        /* number of labels in input data */
  57. int automatic_abscissas=0;
  58. int abscissa_arguments=0;
  59. int curve=0;        /* nonzero if general curve permitted */
  60. int x_arguments=0;
  61. int n;                /* number of entries in x, y, yp, and temp */
  62. int index_array[MAXLABELS];    /* indexes into x and y */
  63. int *p_data=index_array;
  64. int total=100;
  65.  
  66. FILE file;
  67.  
  68. char *label_array[MAXLABELS];
  69. char **p_text=label_array;
  70.  
  71. main(argc,argv) int argc; char **argv;
  72. {    int nn,origin;
  73.     read_data(argc,argv);
  74.     if(breaking && labels)
  75.         {origin=0;
  76.         while(labels--)
  77.             {p_data[0] -= origin;
  78.             nn=p_data[0]+1;
  79.             if(nn) curv0(nn,total);
  80.             origin += nn;
  81.             path += nn;
  82.             x += nn;
  83.             y += nn;
  84.             p_data++;
  85.             p_text++;
  86.             }
  87.         }
  88.     else
  89.         curv0(n,total);
  90. }
  91.  
  92. curv0(n,requested)
  93. int n,requested;
  94. {    int i, j, each, stop, seg=0, regressing=0;
  95.     double s,ds,xx,yy;
  96.  
  97.     for(i=1; i<n; i++) if(path[i]<=path[i-1]) regressing++;
  98.     if (regressing || (curve && (xmin<0. || xmax>1.)))
  99.         {              /* error: echo input */
  100.         for (i=0; i<n; i++)
  101.             {xx=x[i]; if(xlog) xx=exp(xx);
  102.             yy=y[i]; if(ylog) yy=exp(yy);
  103.             printf("%g %g ",xx,yy);
  104.             if(i == p_data[seg]) puts(p_text[seg++]);
  105.             putchar('\n');
  106.             }
  107.         return;
  108.         }
  109.     if(curve) {curv1(path,x,xp,n); curv1(path,y,yp,n);}
  110.     else curv1(x,y,yp,n);
  111.     s=path[0];
  112.     seg=0;
  113.     if(requested<n) requested=n;
  114.     if(x_arguments==2)        /* specific range was requested */
  115.         {s=xmin;
  116.         ds=(xmax-xmin)/requested;
  117.         if(curve) {ds*=path[n-1]; s*=path[n-1];}
  118.         for(i=requested+1; i; i--)
  119.             {value(s,&xx,&yy,n);
  120.             printf("%g %g ",xx,yy);
  121.             if(j==p_data[seg]) puts(p_text[seg++]);
  122.             putchar('\n');
  123.             s+=ds;
  124.             }
  125.         }
  126.     else            /* spline for entire curve was requested */
  127.         {for (i=j=0; j<n-1; j++)
  128.             {stop=requested*(j+1)/(n-1);
  129.             ds=(path[j+1]-path[j])/(stop-i);
  130.             for ( ; i<stop; i++)
  131.                 {value(s,&xx,&yy,n);
  132.                 printf("%g %g ",xx,yy);
  133.                 if(j==p_data[seg]) puts(p_text[seg++]);
  134.                 putchar('\n');
  135.                 s+=ds;
  136.                 }
  137.             }
  138.             xx=x[n-1]; if(xlog) xx=exp(xx);
  139.             yy=y[n-1]; if(ylog) yy=exp(yy);
  140.             printf("%g %g ",xx,yy);
  141.             if(j==p_data[seg]) puts(p_text[seg++]);
  142.             putchar('\n');
  143.         }
  144. }
  145.  
  146. value(s,q,r,n) double s,*q,*r; int n;
  147. {    if(curve)
  148.         {*q=curv2(path,x,xp,s,n);
  149.         *r=curv2(path,y,yp,s,n);
  150.         }
  151.     else
  152.         {*q=s;
  153.         *r=curv2(x,y,yp,s,n);
  154.         }
  155.     if(xlog) *q=exp(*q);
  156.     if(ylog) *r=exp(*r);
  157. }
  158.  
  159. read_data(argc,argv) int argc; char **argv;
  160. {    int i,j,length;
  161.     double xx,yy,d,*pd,sum;
  162.     char *s,*t;
  163. #define BUFSIZE 200
  164.     static char buf[BUFSIZE];
  165.  
  166.     x=path=malloc(ENTRIES*sizeof(double));
  167.     y=malloc(ENTRIES*sizeof(double));
  168.     if(x==0 || y==0) {fprintf(stderr,"can\'t allocate buffer"); exit();}
  169.     if(argc>1 && streq(argv[1],"?")) help();
  170.     if(argc<=1 || *argv[1]=='-') file=stdin;
  171.     else
  172.         {if(argc>1)
  173.             {file=fopen(argv[1],"r");
  174.             if(file==0) {printf("file %s not found\n",argv[1]); exit();}
  175.             argc--; argv++;
  176.             }
  177.         else help();
  178.         }
  179.     argc--; argv++;
  180.     while(argc>0)
  181.         {i=get_parameter(argc,argv);
  182.         argv=argv+i; argc=argc-i;
  183.         }
  184.     if(xlog && !curve)
  185.         {if(x_arguments>1) xmax=mylog(xmax);
  186.         if(x_arguments>=1) xmin=mylog(xmin);
  187.         }
  188.     if(automatic_abscissas && abscissa_arguments<2 && x_arguments>=1)
  189.         abscissa=xmin;
  190.     p_data[0]=-1;
  191.     i=0;
  192.     while(i<ENTRIES)
  193.         {if(fgets(buf,BUFSIZE,file)==0) break;
  194.         t=buf;
  195.         while(*t && isspace(*t)) t++;
  196.         if(*t == 0) continue;        /* skip blank lines */
  197.         buf[strlen(buf)-1]=0;        /* zero out the line feed */
  198.         if(buf[0]==';') {printf("%s\n",buf); continue;}  /* skip comment */
  199.         if(automatic_abscissas)
  200.             {x[i]=abscissa;
  201.             abscissa+=abscissa_step;
  202.             sscanf(buf,"%F",&y[i]);
  203.             if(ylog) y[i]=mylog(y[i]);
  204.             }
  205.         else
  206.             {sscanf(buf,"%F %F",&x[i],&y[i]);
  207.             if(xlog) x[i]=mylog(x[i]);
  208.             if(ylog) y[i]=mylog(y[i]);
  209.             }
  210.         s=buf;                      /* start looking for label */
  211.         while(*s==' ')s++;        /* skip first number */
  212.         while(*s && (*s!=' '))s++;
  213.         if(!automatic_abscissas)    /* skip second number */
  214.             {while(*s==' ')s++;
  215.             while(*s && (*s!=' '))s++;
  216.             }
  217.         while(*s==' ')s++;
  218.         if((length=strlen(s))&&(labels<MAXLABELS))
  219.             {if(*s=='\"')           /* label in quotes */
  220.                 {t=s+1;
  221.                 while(*t && (*t!='\"')) t++;
  222.                 t++;
  223.                 }
  224.             else                    /* label without quotes */
  225.                 {t=s;
  226.                 while(*t && (*t!=' '))t++;
  227.                 }
  228.             *t=0;
  229.             length=t-s;
  230.             p_data[labels]=i;
  231.             p_text[labels]=malloc(length+1);
  232.             if(p_text[labels]) strcpy(p_text[labels++],s);
  233.             }
  234.         i++;
  235.         }
  236.     if(breaking && (!labels || p_data[labels-1]!=n-1))
  237.         {p_data[labels]=i-1;
  238.         if(p_text[labels]=malloc(1)) *p_text[labels++]=0;
  239.         }
  240.     n=i;
  241.     yp=malloc(n*sizeof(double));
  242.     temp=malloc(n*sizeof(double));
  243.     if(temp==0 || yp==0) {fprintf(stderr,"can\'t allocate buffer"); exit();}
  244.     if(curve)
  245.         {xp=malloc(n*sizeof(double));
  246.         path=malloc(n*sizeof(double));
  247.         if(xp==0|| path==0) {fprintf(stderr,"can\'t allocate buffer"); exit();}
  248.         path[0]=sum=0.;
  249.         for (i=1; i<n; i++)
  250.             {xx=x[i]-x[i-1];
  251.             yy=y[i]-y[i-1];
  252.             path[i]=(sum+=sqrt(xx*xx + yy*yy));
  253.             }
  254. /*        for(i=0; i<n; i++)
  255.             printf("path[%d]=%g  x[%d]=%g \n",i,path[i],i,x[i]); */
  256.         }
  257. }
  258.  
  259.  
  260. /* get_parameter - process one command line option
  261.         (return # parameters used) */
  262. get_parameter(argc,argv) int argc; char **argv;
  263. {    int i;
  264.     if(streq(*argv,"-a"))
  265.         {i=get_double(argc,argv,2,&abscissa_step,&abscissa,&abscissa);
  266.         abscissa_arguments=i-1;
  267.         automatic_abscissas=1;
  268.         return i;
  269.         }
  270.     else if(streq(*argv,"-b")) {breaking=1; return 1;}
  271.     else if(streq(*argv,"-c")) {curve=1; return 1;}
  272.     else if(streq(*argv,"-d")) {debugging=1; return 1;}
  273.     else if(streq(*argv,"-n"))
  274.         {if((argc>1) && numeric(argv[1])) total=atoi(argv[1]);
  275.         if(debugging)printf("%d output pairs",total);
  276.         return 2;
  277.         }
  278.     if(streq(*argv,"-t"))
  279.         {return(get_double(argc,argv,1,&sigma,&abscissa,&abscissa));
  280.         }
  281.     else if(streq(*argv,"-x"))
  282.         {i=get_double(argc,argv,2,&xmin,&xmax,&xmax);
  283.         x_arguments=i-1;
  284.         return i;
  285.         }
  286.     else if(streq(*argv,"-xl")) {xlog++; return 1;}
  287.     else if(streq(*argv,"-yl")) {ylog++; return 1;}
  288.     else gripe(argv);
  289. }
  290.  
  291. get_double(argc,argv,permitted,a,b,c)
  292. int argc,permitted; char **argv; double *a,*b,*c;
  293. {    int i=1;
  294.     if((permitted--)>0 && (argc>i) && numeric(argv[i])) *a=atof(argv[i++]);
  295.     if((permitted--)>0 && (argc>i) && numeric(argv[i])) *b=atof(argv[i++]);
  296.     if((permitted--)>0 && (argc>i) && numeric(argv[i])) *c=atof(argv[i++]);
  297.     return i;
  298. }
  299.  
  300. int streq(a,b) char *a,*b;
  301. {    while(*a)
  302.         {if(*a!=*b)return 0;
  303.         a++; b++;
  304.         }
  305.     return 1;
  306. }
  307.  
  308. gripe_arg(s) char *s;
  309. {    fprintf(stderr,"argument missing for switch %s",s);
  310.     help();
  311. }
  312.  
  313. gripe(argv) char **argv;
  314. {    fprintf(stderr,*argv); fprintf(stderr," isn\'t a legal argument \n\n");
  315.     help();
  316. }
  317.  
  318. help()
  319. {    fprintf(stderr,"spline   version %s",VERSION);
  320.     fprintf(stderr,"\nusage: spline [file][options]\n");
  321.     fprintf(stderr,"options are:\n");
  322.     fprintf(stderr,"  -a  [step [start]] automatic abscissas \n");
  323.     fprintf(stderr,"  -b                 break interpolation after each label \n");
  324.     fprintf(stderr,"  -c                 general curve \n");
  325.     fprintf(stderr,"  -n  num            interpolate over num intervals \n");
  326.     fprintf(stderr,"  -t  num            tension in interpolating curve \n");
  327.     fprintf(stderr,"  -x  [min [max]]    interpolate from min to max \n");
  328.     fprintf(stderr,"  -xl                take logs of x values before interpolating \n");
  329.     fprintf(stderr,"  -yl                take logs of y values before interpolating \n");
  330.     exit();
  331. }
  332.  
  333. numeric(s) char *s;
  334. {    char c;
  335.     while(c=*s++)
  336.         {if((c<='9' && c>='0') || c=='+' || c=='-' || c=='.') continue;
  337.         return 0;
  338.         }
  339.     return 1;
  340. }
  341.  
  342. curv1(x,y,yp,n) double *x,*y,*yp; int n;
  343. {    int i;
  344.     double c1,c2,c3,deln,delnm1,delnn,dels,delx1,delx2,delx12;
  345.     double diag1,diag2,diagin,dx1,dx2,exps;
  346.     double sigmap,sinhs,sinhin,slpp1,slppn,spdiag;
  347.  
  348.     delx1=x[1] - x[0];
  349.     dx1=(y[1] - y[0])/delx1;
  350.     if(sigma >= 0.) {slpp1=slp1; slppn=slpn;}
  351.     else
  352.         {if(n!=2)
  353.             {delx2= x[2] - x[1];
  354.             delx12= x[2] - x[0];
  355.             c1= -(delx12 + delx1)/delx12/delx1;
  356.             c2= delx12/delx1/delx2;
  357.             c3= -delx1/delx12/delx2;
  358.             slpp1= c1*y[0] + c2*y[1] + c3*y[2];
  359.             deln= x[n-1] - x[n-2];
  360.             delnm1= x[n-2] - x[n-3];
  361.             delnn= x[n-1] - x[n-3];
  362.             c1= (delnn + deln)/delnn/deln;
  363.             c2= -delnn/deln/delnm1;
  364.             c3= deln/delnn/delnm1;
  365.             slppn= c3*y[n-3] + c2*y[n-2] + c1*y[n-1];
  366.             }
  367.         else yp[0]=yp[1]=0.;
  368.         }
  369.             /* denormalize tension factor */
  370.     sigmap=fabs(sigma)*(n-1)/(x[n-1]-x[0]);
  371.             /* set up right hand side and tridiagonal system for
  372.                yp and perform forward elimination    */
  373.     dels=sigmap*delx1;
  374.     exps=exp(dels);
  375.     sinhs=.5*(exps-1./exps);
  376.     sinhin=1./(delx1*sinhs);
  377.     diag1=sinhin*(dels*.5*(exps+1./exps)-sinhs);
  378.     diagin=1./diag1;
  379.     yp[0]=diagin*(dx1-slpp1);
  380.     spdiag=sinhin*(sinhs-dels);
  381.     temp[0]=diagin*spdiag;
  382.     if(n!=2)
  383.         {for(i=1; i<=n-2; i++)
  384.             {delx2=x[i+1]-x[i];
  385.             dx2=(y[i+1]-y[i])/delx2;
  386.             dels=sigmap*delx2;
  387.             exps=exp(dels);
  388.             sinhs=.5*(exps-1./exps);
  389.             sinhin=1./(delx2*sinhs);
  390.             diag2=sinhin*(dels*(.5*(exps+1./exps))-sinhs);
  391.             diagin=1./(diag1+diag2-spdiag*temp[i-1]);
  392.             yp[i]=diagin*(dx2-dx1-spdiag*yp[i-1]);
  393.             spdiag=sinhin*(sinhs-dels);
  394.             temp[i]=diagin*spdiag;
  395.             dx1=dx2;
  396.             diag1=diag2;
  397.             }
  398.         }
  399.     diagin=1./(diag1-spdiag*temp[n-2]);
  400.     yp[n-1]=diagin*(slppn-dx2-spdiag*yp[n-2]);
  401.                     /* perform back substitution */
  402.     for (i=n-2; i>=0; i--) yp[i] -= temp[i]*yp[i+1];
  403. }
  404.  
  405.  
  406. double curv2(x,y,yp,t,n) double *x,*y,*yp,t; int n;
  407. {    int i,j;
  408.     static int i1=1;
  409.     double del1,del2,dels,exps,exps1,s,sigmap,sinhd1,sinhd2,sinhs;
  410.  
  411.     s=x[n-1]-x[0];
  412.     sigmap=fabs(sigma)*(n-1)/s;
  413. #ifdef WORK
  414.     for (j=2; j; j--)    /* want: x[i-1] <= t < x[i], 0 < i <= n */
  415.         {for (i=i1; i<n; i++)
  416.             {if(x[i]>t) break;
  417.             }
  418.         if(i==n) i=n-1;
  419.         if(x[i-1]<=t || t<=x[0]) break;
  420.         i1=1;
  421.         }
  422. #endif
  423.     i=i1;
  424.     while(i<n && t>=x[i]) i++;
  425.     while(i>1 && x[i-1]>t) i--;
  426.     i1=i;
  427.     del1=t-x[i-1];
  428.     del2=x[i]-t;
  429.     dels=x[i] - x[i-1];
  430.     exps1=exp(sigmap*del1); sinhd1=.5*(exps1-1./exps1);
  431.     exps= exp(sigmap*del2); sinhd2=.5*(exps-1./exps);
  432.     exps=exps1*exps;        sinhs=.5*(exps-1./exps);
  433.     return ((yp[i]*sinhd1 + yp[i-1]*sinhd2)/sinhs +
  434.             ((y[i] - yp[i])*del1 + (y[i-1] - yp[i-1])*del2)/dels);
  435. }
  436.  
  437. double mylog(x) double x;
  438. {    if(x>0.) return log(x);
  439.     fprintf(stderr,"can%'t take log of nonpositive number");
  440.     exit(1);
  441.     return 0.;
  442. }
  443.  
  444.